%input: immagine di input
%q_dim: dimensione dei quadrati ce inscrivono i punti
%level_precision: numero di cerchi di raggio diffrente
function output = halftoning_function(input, q_dim, level_precision)

if(level_precision > 256)
    display('level_number automatically set to 256');
    level_precision = 256;
end

if(mod(q_dim, 2) == 0)
    q_dim = q_dim+1;
end

fattore_scala = 5;
input = imresize(input, fattore_scala);

sovrapposizione = ceil(q_dim*(1 - sqrt(3)/2)-1);

[h,w] = size(input);
output = ones(h, w)*255;

zircle = imread('dot64.png');
zircle = imresize(zircle, [q_dim, q_dim]);


posx = -(floor(q_dim/2));
posy = -(floor(q_dim/2));
line_counter = 0;

while (posy <= h)
    line_counter = line_counter+1;
    if(mod(line_counter, 2) == 0)
        posx = 1;
    else
        posx = -(floor(q_dim/2));
    end
    
    while(posx <= w)
        %devo controllare se il quadrato esce o meno dai confini dell
        %immagine
        inizio_up = 1;
        fine_down = q_dim;
        inizio_left = 1;
        fine_right = q_dim;
        
        %controllo se esce sopra l'immagine
        if(posy < 1)
            inizio_up = 1-posy;
        end
        %controllo se esce sotto l'immagine
        if((posy + fine_down) > h)
            fine_down = h - posy;
        end
        
        %controllo se esce da sinistra
        if(posx < 1)
            inizio_left = 1-posx;
        end
        %controllo se esce a destra dell'immagine
        if((posx + fine_right) > w)
            fine_right = w - posx;
        end
        
        average = 0;
        for i = inizio_up : fine_down
            for j = inizio_left : fine_right
                %if(output(posy+i, posx+j) == 255)
                    %output(posy+i, posx+j) = zircle(i, j);
                %end
                average = average + double(input(posy+i, posx+j));
            end
        end
        
        average = round(average/(q_dim*q_dim));
        dimension = circle_dimension(level_precision, q_dim, average);
        cerchio = circle_of_dimension(zircle, dimension, q_dim);
        for i = inizio_up : fine_down
            for j = inizio_left : fine_right
                if(output(posy+i, posx+j) == 255)
                    output(posy+i, posx+j) = cerchio(i, j);
                end
            end
        end
        posx = posx + q_dim;
    end
    posy = posy + q_dim - sovrapposizione;
end

%output = imresize(output, 1/fattore_scala);
[h, w] = size(output);
for i = h
    for j = w
        if(output(i,j) < 250)
            output(i,j) = 0;
        else
            output(i,j) = 255;
        end
    end
end

end

%ritorna la dimensione che deve avere il cerchio in base al livello
%di precisione e alla media
function dimension = circle_dimension(level_precision, q_dim, average)
%x = (level_precision*floor(average/level_precision));
%y = x*(-q_dim/255) + q_dim;
%y = round(level_precision*(y/level_precision));
dimension = round(average*(-q_dim/255) + q_dim);
if(dimension < 0)
    dimension = 0;
elseif(dimension > q_dim)
    dimension = q_dim;
end
if(mod(dimension, 2) == 0 && dimension ~= 0)
    dimension = dimension+1;
end
end

%restituisce un cerchio nero di diametro dimension inscritto
%in un quadrato di dimensione q_dim a sfondo binco
function cerchio = circle_of_dimension(circle, dimension, q_dim)
if(dimension == q_dim)
    cerchio = imresize(circle, [q_dim, q_dim]);
elseif(dimension == 0)
    cerchio = ones(q_dim, q_dim)*255;
else
    if(mod(dimension, 2) == 0)
        dimension = dimension+1;        
    end
    cerchio = imresize(circle, [dimension dimension]);
    A = uint8(ones(floor((q_dim-dimension)/2), dimension))*255;
    B = uint8(ones(q_dim, floor(q_dim - dimension)/2)) * 255;
    cerchio = vertcat(A, cerchio, A);
    cerchio = horzcat(B, cerchio, B);
end
end